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Abstract 

Numerical micropermeametry is performed on three dimensional porous samples having a linear 
size of approximately 3 mm and a resolution of 7.5 urn. One of the samples is a microtomographic 
image of Fontainebleau sandstone. Two of the samples are stochastic reconstructions with the same 
porosity, specific surface area, and two-point correlation function as the Fontainebleau sample. The 
fourth sample is a physical model which mimics the processes of sedimentation, compaction and 
diagenesis of Fontainebleau sandstone. The permeabilities of these samples are determined by nu- 
merically solving at low Reynolds numbers the appropriate Stokes equations in the pore spaces of 
the samples. The physical diagenesis model appears to reproduce the permeability of the real sand- 
stone sample quite accurately, while the permeabilities of the stochastic reconstructions deviate 
from the latter by at least an order of magnitude. This finding confirms earlier qualitative pre- 
dictions based on local porosity theory. Two numerical algorithms were used in these simulations. 
One is based on the lattice-Boltzmann method, and the other on conventional finite-difference 
techniques. The accuracy of these two methods is discussed and compared, also with experiment. 

PACS: 81.05.Rm, 83.85.Pt, 61.43.Gt 
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I. INTRODUCTION 



Almost all investigations of porous media focus on the prediction of effective material 
properties such as fluid permeability, electric or thermal conductivity, or elastic constants 

|2|. The knowledge or at least a reliable prediction of these properties is of great interest 
in a wide field of technical applications ranging from petroleum engineering || |J , to paper 
manufacturing and contaminant transport ||. The predictions are obtained either from 
approximate theories which link the physical properties to geometrical observables, from 
geometrical models for which the physical problem can be solved more easily, or from various 
cross property relations, which relate the parameter in question to other physical transport 
parameters. 

In this context, the exact numerical calculation of transport properties serves three pur- 
poses: (i) testing and validation of theories and theoretical predictions, (ii) comparison of 
geometrical models, and (iii) testing of faithfulness of computerized tomographic imaging 
by comparing numerically calculated transport parameters with their experimental values. 

Of particular interest for porous media is the permeability, and more precisely its fluctua- 
tions. These fluctuations are important because they dominate the large scale permeability. 
For this reason it is important to collect as many micropermeametry measurements as pos- 
sible. Experimental micropermeametry is costly and inaccurate. Hence exact numerical 
calculations are becoming an interesting alternative for studying fluctuations in permeabil- 
ity. 

In this article we compare the permeabilities of a three-dimensional computerized tomo- 
graphic image of Fontainebleau sandstone and its three physical and stochastic reconstruc- 
tion models. We find that especially the stochastic models fail in reconstructing the fluid 
permeability of the original sandstone. This finding is in good agreement with conclusions 
drawn previously from a purely geometrical characterization of the same microstructures 

0- 

The exact numerical calculation of transport parameters for digitized three-dimensional 
samples remains a computationally demanding task, and only few studies exist which test 
the accuracy of such calculations. Here, we compare the results obtained by means of a 
finite difference (FD) method and a lattice-Boltzmann (LB) algorithm. We begin our study 
by calibrating both simulation methods against exact solutions of the Stokes equation for 
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straight tubes and cubic arrays of spheres. These calculations also serve to compare the 
speed of both methods. We then proceed to apply both algorithms to the experimental 
sample and its three models. 



II. DEFINITION OF THE PROBLEM 

The problem to be solved is that of slow laminar flow through a three-dimensional porous 
medium on a microscopic level. The three-dimensional microstructure of a two-phase porous 
medium § consisting of a pore phase P and a matrix or rock phase M with S = P U M is 
described in detail by the characteristic function \ of a single phase G G {M, P} with 

{1 for feG 
(1) 
for x £ G 

In the following x is the position vector of a cubic lattice, x = a ■ X\e\ + a ■ X2&2 + a ■ x^, 
with Xi = 0, 1, . . . , Mi — 1, the unit vectors e*j of the Cartesian coordinate system, and the 
grid spacing a. The total number of lattice points is given by N = M1M2M3. 

The Reynolds numbers of interest in geophysical and petrophysical applications are usu- 
ally much smaller than unity || and hence it suffices to solve the Stokes equation. In the 
pore space geometry described by the characteristic function x p (%)i the steady-state Stokes 
equation and the constraint of incompressibility read 

T]Av(x) - Vp(f) = 0; x G P, (2) 
V ■ v(x) =0; ?GP . (3) 

On the pore-matrix interface <9P we apply no-flow boundary conditions, 

v(x) = 0; x G d P . (4) 

Darcy's law permits us to compute the j-th column of the macroscopic permeability tensor 
from the microscopic solution of the hydrodynamic problem Eqs. (0) - @) for a pressure 
gradient applied along the direction according to 

t v = (vm^^^ ■ (5) 

Pin Pont 

Here the pressures at the inlet and the outlet surface are given by p m and p out as defined 
below in Eqs. flTBD and (0) |4l| , respectively, (. . .)^ g § denotes an average over all lattice 



points, and Vi = v ■ e*j. 



III. SAMPLES 



We investigate four different samples §[=X' ^DM ' ^GF' an< ^ *%A- ^ ne ^ rs ^ sam pl e §>EX> or 
in abbreviated form EX, was obtained experimentally by means of computerized tomography 
from a core of Fontainebleau sandstone. This sandstone is a popular reference standard 
because of its exceptional chemical, crystallographic and microstructural simplicity 10 



Fontainebleau sandstone consists of crystalline quartz grains that have been eroded for long 
periods before being deposited in dunes along shore lines during the Oligocene, i.e. roughly 
30 million years ago. It is well sorted containing grains of around 200 fim in diameter. During 
its geological evolution, which is still not fully understood, the sand was cemented by silica 
crystallizing around the grains. Fontainebleau sandstone exhibits intergranular porosity 
ranging from 0.03 to roughly 0.3 [fLOl - The computer-assisted microtomography was carried 
out on a micro-plug drilled from a larger original core. The original core from which the 
micro-plug was taken had porosity <fi* = 0.1484, permeability k* = 1.3 D, (1 D=0.987 /im 2 ) 
and formation factor 22.1 (dimensionless electrical resistivity [|ll|]). The microtomographic 
dataset has dimension Mi x M 2 x M 3 = 299 x 300 x 300 with a resolution of a = 7.5 /im, 
and porosity <fi = 0.1355. The pore space P^y is visualized in Fig. 1 of Ref. 0. 

The three remaining samples are physical and stochastic reconstruction models for the 
Fontainebleau sample EX. All have the same lattice resolution, a = 7.5/xm, and approxi- 
mately the same porosity. The "diagenesis model" DM tries to mimic the geological forma- 
tion process of the natural sandstone in three steps: the sedimentation of spherical grains, 
the compaction of the sediment, and the simulation of quartz cement overgrowth. The 
sample dimensions are Mi x M 2 x M 3 = 255 x 255 x 255. 

The SA and GF samples are stochastic models with dimensions Mi x M 2 x M3 = 256 x 
256 x 256. Both models reconstruct the porosity <fi and the two-point correlation function 
of the original sandstone EX. This implies the reconstruction of the specific surface Sy. 
However, due to problems in the reconstruction procedure of the Gaussian field method, the 
porosity and the specific surface of the GF model do not match exactly those of the original 
sandstone. We find = 0.1354 for SA, and = 0.1421 for GF. 

For a more detailed description of the modeling procedures and a visualization of the 
microstructures, the reader is referred to Ref. J?]] and the references therein. In the same 
article results of an extensive geometrical investigation of the four samples are presented, 
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which use both classical geometric quantities but also concepts introduced in local porosity 
theory |L2], [13]. The main findings in Ref. were as follows: 



1. None of the three models can reproduce the visual appearance of the original sandstone. 
EX shows a granular structure of the matrix phase where single sand grains can easily 
be identified. The matrix phase of DM is also clearly granular but with artificial, 
spherical grains. In both models the matrix and the pore phase are very well connected, 
and the pore-matrix interface is smooth. In contrast, the pore-matrix interface of the 
stochastic models SA and GF is very rough. Here, both phases are strongly scattered 
and exhibit isolated clusters. 

2. The two-point correlation functions, S2<r) = (x(^i)x(^2)) with r = \x% — X2I, of SA 
and GF show good agreement with the original sandstone except for minor deviations 
at small r in the case of GF. The correlation function of DM clearly deviates from that 
of EX. Moreover, it shows strong anisotropy with respect to directions e*i, e*2 and e^. A 
modified correlation function giving the conditional probability to find two points in 
pore space that are also connected by a path inside the pore space was measured and 
discussed in Ref. [[HJ . Small but significant differences exist between the samples. The 
experimental sample EX is more stable under the morphological operations of erosion 
and dilation |T5], |TjJ. 

3. The differences between the samples are most pronounced when comparing the geo- 
metrical connectivity of the pore space. As a measure for the geometrical connectivity, 
we use the total fraction of percolating cells Pz{L) at scale L, introduced in local poros- 
ity theory which is defined as the probability for a cubic subblock of size L of the 
sample to percolate in all three coordinate directions e*j. Here, percolation in direction 
e*j means that there exists a path lying entirely in the pore space, which spans from 
one face of the cubic subsample perpendicular to e*j to the opposite face. For EX and 
DM the curves of Ps{L) nearly coincide while for the stochastic models SA and GF the 
curves of Ps(L) fall well below that of the original sandstone (see Fig. 13 in Ref. 0). 
Table | gives the values of Ps{L) for L = 60a. Again we find anisotropy in the DM 
model when we measure the probability pg i (L) for a subblock to percolate in directions 



Geometrical connectivity is an indispensable precondition for dynamical connectivity and 
physical transport. Hence, we expect to find a strong correlation between the total fraction 
of percolating cells p$(L) and the macroscopic permeability whose calculation is discussed 
next. 



IV. THE NUMERICAL METHODS 



A. The finite difference method 



1. Algorithm 

Numerically we obtain the solution of Eqs. (0) - (01) from the infinite time limit of the 



time dependent Stokes problem using an iterative pressure-correction algorithm fl4| , |17 . 
Discretization in time of the time-dependent Stokes equation yields 

V"+\X)-^(X) = _ V^+l(f), ( 6 ) 

V ■ iT +1 (f) = 0, (7) 

where the superscript n denotes the iteration step. In our case the discretized time derivative 
on the left hand side of Eq. (|6|) has no physical meaning. In the long-time limit t7 ra+1 (af) = 
v n {x) holds, and we recover Eq. (§). 

Given the solutions v n and p n at iteration step n, an approximate solution v* for the 
velocity field is obtained from 



At 

Subtracting Eq. (H) from Eq. (H), we find that 



r/AiT l (f) - Vp n (x) . (8) 



v«+\x)-v*{x) _ ^/ n+1 ^ 



A/ -V(p n+ \x)-p n {x)) . (9) 

On the right hand side of this equation appears the pressure correction p'{x) = p n+1 (x) — 
p n (x). Applying the V operator to Eq. @, and using the incompressibility constraint Eq. 
(0) we obtain Poisson's equation for p', 

Ap'(x) = — V • v* (x) . (10) 

Thus, we arrive at the following algorithm: 



1. Let v" 1 and p n be the solution of the velocity and the pressure field, respectively, at 
iteration step n with the maximum absolute error e n = max^p [rjAv 11 ^) — Vp(x)|. 
From v" - and p n an approximate solution v* of the velocity field is calculated using Eq. 

®- 

2. Using the definition of the pressure correction, a new pressure p n+l = p n + p' is 
obtained from a solution of Eq. fllPp. This part of the algorithm consumes most of the 
computation time, because Eq. ( |10D has to be solved for each iteration step. However, 
we found that it suffices to solve Eq. ((TOj) only up to an error 

max \Ap'(x) — — V • v*(x)\ < ^e n 
x& At 

where 7 is an empirical factor. For the calculations presented here, a value of 7 in the 
range 0.01 < 7 < 0.1 seems to be appropriate. 

We use a successive over-relaxation method to solve Eq. (]10|) . Of course it would be 
desirable to use more sophisticated methods such as e.g. a multigrid method, but we 
could not find a general procedure to restrict the microstructure of the porous medium 
to a coarser grid without changing the topology of the pore space. 

3. From p' and v* a new velocity v"- +1 is calculated using Eq. (|S]). The algorithm 
terminates when e n+1 is smaller than some given value. 



The equations are spatially discretized using a marker-and-cell (MAC) grid |18fl . The 
pressure values are placed at the centers of the grid cells. They coincide with the lattice 
points of the discretized characteristic function x v - On each face of a grid cell the velocity 
component perpendicular to this face is located. The pore-matrix interface dP follows the 
surface of the cubic grid cells. For velocity components perpendicular to the interface the 
boundary condition Eq. (^), 

v±(x) = 0, (11) 

is implemented exactly. For parallel velocity components at distance a/2 from the interface, 
Eq. (f|) is fulfilled to second order accuracy, 

V \\(x) = -v\\(x + ae ± ) + 0(a 2 ), (12) 
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in which the interface is located as shown in Fig. |I[ Inserting Eq. (|j) into the Stokes 
equation Eq. (H) one finds as the boundary conditions for the pressure 



d±p{x)\ dP = . 

Using Eqs. ([□]) - (|13"|), the A operators in Eqs. (|]) and (|T0| ) take the form 
Avi(x) — Vi(x — ati) + Vi(x + ati) — 2vi(x) 

+ ^2 x f^ + a ^')^p(^ + + [vi(x + aej) + Vi(x)] 

+ X P (x ~ aej)x P (£ - ae 3 - + aci) [vi(x - aej) + Vi(x)] - Avi(x) 



x e P, 



and 



Ap'(x) = 



X p (x + ae { ) (p'(x + aii) - p'(x)) 
X p (x - ae % ) (p'(x - aei) - p'(x)) 



■x e P, 



(13) 



(14) 



(15) 



where we used second order accurate central differences to discretize the spatial derivatives. 



2. Boundary conditions 

On the sample surface an additional outer layer of grid cells, the so-called shadow row, 
is added. It provides the neighboring pressure and velocity values needed for the evaluation 
of Eqs. (|l"4"D and (|15D for the grid cells on the sample surface. The pressure and velocity 
values of the shadow row are set according to the macroscopic boundary conditions. Let e*,- 
denote the direction of the applied pressure gradient. For grid cells of the shadow row on 
the outflow boundary of the sample, i.e. x G {x : Xj = Mj}, we choose 

p(x) = p ou t, (16) 
Vi(x) =0 for i 7^ j, (17) 
Vj(x) = Vj(x — ae\\), (18) 

as the boundary conditions. For x G {x : Xj = —1}, i.e. on the inflow boundary, we set 

P(x) = Pin, (19) 



Pin, 

Vi(x) =0 for i ^ j, 
v j(x) = Vj(x — ae»), 



(20) 
(21) 
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as the boundary conditions. The last condition has to be introduced directly into Eq. fll4"D 
because the position x — aej lies outside the simulation lattice. In the simulations reported 
below we used always p in = 1, p out = — 1 unless indicated otherwise. 

On the remaining sample surfaces the grid cells of the shadow row are assigned to the 
matrix phase. The velocity and pressure values are set to zero. 



B. The lattice-Boltzmann method 



In this section we introduce the lattice-Boltzmann method used here, in particular the 
LBGK (lattice-Bhatnagar-Gross-Krook) model. We discuss then the basic hydrodynamics 
of the model and the relevant boundary conditions. The numerical accuracy of the lattice- 
Boltzmann results for permeability will be considered in terms of finite-size effects. 



1. Lattice-Boltzmann hydrodynamics 



The lattice-Boltzmann method [19, ^(J, is a mesoscopic approach for computational 



fluid dynamics in which the basic idea is to solve a discretized Boltzmann equation. The 
macroscopic dynamics of the system can be shown to obey the Navier-Stokes equation. One 
of the most successful applications of the method has been to flow in porous media [[5], ^3|. 24 



In this method the fluid is modeled by particle distributions that move on a regular 
lattice. In our implementation each lattice point is connected to its nearest and next nearest 
neighbors. Together with a rest particle, each lattice point is then occupied by 19 different 
particles (the D3Q19 model). At each time step particles propagate to their adjacent lattice 
points, and re-distribute their momenta in the subsequent collisions. The dynamics of the 
LBGK model is given by the equation ITDL 120 



/<(r + c,,f + 1) = fi(v,t) + ~[jT(r,i) - (22) 

where c, is a vector pointing to an adjacent lattice site, fi(r, t) is the density of the particles 
moving in the c r direction, r is the BGK relaxation parameter, and / 4 eq (r, t) is the equilibrium 
distribution towards which the particle populations are relaxed. Hydrodynamic quantities 
like density p and velocity u are obtained from the velocity moments of the distribution /j 
in analogy with the kinetic theory of gases. The equilibrium distribution can be chosen in 
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many ways. A common choice is 

/r = U(l + ^(c, • u) + ^(c t • u) 2 - ±- 2 u\ (23) 

in which tj is a weight factor that depends on the length of the link vector Cj, and c s is the 
speed of sound in the fluid. The weights t« we choose here in accordance with the 19-link 
LBGK model, and they are ~, and ^ for the rest particle and the particles moving to 
the nearest and next-nearest neighbor sites, respectively. The speed of sound is c s = ^= for 
this model, and the kinematic viscosity of the simulated fluid is rj = 2j -^-- (Here and in the 
following, lattice units are always used if the units are not specified.) The fluid pressure is 
given by 

p(r, t) = c 2 (p(r, t)-p) = c 2 Ap(r, t), (24) 

where p is the mean density of the fluid. 

The Stokes equation, Eq. (^), is produced directly by the linearized lattice-Boltzmann 
method, in which the quadratic velocity terms in the equilibrium distribution function Eq. 
(f23|) are neglected. To be consistent with the finite difference method, we use in what follows 
the linearized lattice-Boltzmann method if not stated otherwise. 



2. Boundary conditions 

The physical boundary condition at solid-fluid interfaces is the no-flow condition Eq. (|]), 
which in lattice-Boltzmann simulations is usually realized by the so-called bounce-back rule 



25| , p26f . In this approach the momenta of the particles that meet a solid wall are simply 
reversed. 

In simple shear flows the bounce-back condition assumes that the location of the wall is 
exactly halfway between the last fluid point and the first wall point. In more complicated 
cases the no-flow boundary lies somewhere in between these two points, the exact place 
depending on the relaxation parameter and the geometry of the system |26| , |27|j . In Poiseuille 
flow, e.g., the bounce-back rule gives velocity fields that deviate from the exact solution, for 
no-flow boundaries at exactly halfway between the last fluid point and the first solid point, 

by m 



A - 48r/ 2 -477-I 

LAU — lt s im ^exact — ^max ' 
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where w S i m and M exa ct are the simulated and the exact velocities, respectively, M max is the 
velocity at the center of the channel, and L is channel width. This implies that the simu- 
lated permeability will depend somewhat on viscosity especially at low discretization levels. 
This viscosity dependence can practically be eliminated by using the so-called second or- 
der boundaries, in which case the desired location of the no-flow boundary is determined by 
extrapolating the distribution function from the last fluid points. Some of these more sophis- 



ticated solid- fluid boundaries are restricted to regular geometries |28, 29j, but there are also 



general boundary-fitted models p0| , pl|] available. For practical simulations the bounce-back 
boundary is however very attractive, because it is a simple and computationally efficient 
method for imposing no-flow conditions on irregularly shaped walls. Also, the error created 
by the bounce-back boundary does not destroy the spatial second-order convergence of the 



method ^7, |2 



In simulating fluid flow it is important that the velocity and pressure boundary conditions 
of the system have been imposed in a consistent way. However, general velocity and pressure 



boundaries are still under development for the lattice-Boltzmann method |33L |34| . So far 



in most of the practical simulations a body force has been implemented |23|, |2^, [35| instead 
of pressure or velocity boundaries. 

When the body force is used, the pressure gradient acting on the fluid is replaced with 
a uniform external force. The use of a body force is based on the assumption that, on 
the average, the effect of an external pressure gradient is constant throughout the system, 
and that it can thus be replaced with a constant force that adds at every time step a fixed 
amount of momentum on the fluid points. Conditions that are close to pressure boundaries 
can be obtained by averaging the velocity and pressure fields over the planes of the inlet 
and outlet of the simulated system |27| . 



During one iteration step, the fluid momentum oscillates in the stationary state by an 
amount given to each fluid point by the body force. For this reason the fluid velocity is now 



defined as the average of the pre-collision and post-collision values []36], 37 1 



Pressure fields generated by the body force are obtained from the effective pressure p e g, 

p eff (r, t) = c 2 s Ap(r, t) - pgx, (26) 

where x is the distance from the inlet of the system measured in the flow direction, and g is 
the acceleration the body force gives to the fluid. 
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It is a well-known fact that, due to staggered invariants, the fluid momentum may oscillate 



in a time scale of a few time steps ||22|| , even in the stationary state. In open areas this 
effect is usually unimportant, but in closed pores this effect may become visible as the fluid 
momentum may oscillate around zero, with a magnitude determined by the body force. This 
effect may lead to some corruption of the fluid- velocity distributions as can be seen in Fig|| 
Notice that staggered momenta can be eliminated by averaging the quantities over a few 
time steps. 

Notice finally that the diagonal links allow the fluid to leak to neighboring lattice points 
which have only a single edge in common. For this reason the 'standard' lattice-Boltzmann 
model is not expected to be accurate very close to the percolation threshold. For the 3D 
checkerboard structure, e.g., we found that the permeability of the system was about 0.036 
lattice units for all the six lattice resolutions that were used, although the structure is not 
percolating. If better accuracy is needed, diagonal leaks can be eliminated by applying the 
bounce-back rule on such diagonal links which actually cross a solid boundary (like diagonals 
in the checkerboard structure). 

3. Finite-size effects and the saturation time 

The accuracy of lattice-Boltzmann simulations depends on the ratio of the mean free path 



A m fp of the fluid particles to the representative size Ao of the obstacles and pores |21], [23], |34| . 
The simulated flow field does not describe the true hydrodynamic behavior unless this ratio 
is small. For increasing A m f p / Aq ratio Knudsen-flow behavior is found, which is also true in 



real fluids |38| . These effects must always be considered when lattice-Boltzmann simulations 



are performed. In this way the maximum size of the lattice spacing can be estimated together 
with the accuracy of the simulations. 

Finite-size effects restrict to some extent the use of LB-methods based on regular lattices. 
In porous media close to the percolation threshold, e.g., many pores are very small, and very 
big lattices may be needed for realistic simulations. It is still an open question whether the 
finite-size effects are always dominated by the minimum pore size or the average pore size. 
The effect can be estimated by simulating the system with several different lattice spacings, 
but occasionally it is difficult to distinguish the finite-size effects from other sources of 
numerical error. 
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Practice has shown || ^3] that smaller values of the relaxation parameter r tend to 



decrease the finite-size effects (see also our simulations below). Equation ( p5|) can be used 
to explain this: due to Knudsen-flow effects, low-discretization simulations regularly give 
too high permeabilities, whereas decrease of r has the opposite effect down to r = 0.625, at 



least for tube flows. On the other hand, the lattice-Boltzmann algorithm may become fF2 



unstable with values of r close to 0.5. In practical permeability simulations the relaxation 
parameter has usually been chosen to be bigger than 0.6. The effect of r on the behavior of 
the lattice-Boltzmann model is thus quite complicated and not yet fully understood. 

In permeability simulations a simple dimensional analysis shows that, with a constant 
body force, the saturation time t sat needed to reach the steady state is of the form 

t sat (x R 2 po Jr], (27) 

where -R por e is the characteristic length of the void pores in the system. For systems with 
high porosity, the saturation times can therefore be very long. In some cases, tens of thou- 
sands of time steps may be needed. It is thus evident that a constant body force may be 
computationally inefficient, especially when one is only interested in the steady-state solu- 
tion. The saturation time can be reduced by using e.g. the iterative momentum-relaxation 
(IMR) method, where the applied body force is adjusted during the iteration in a definite 



relation to the change of the fluid momentum during iteration steps [27]. For other ways to 



reduce the saturation time see 39 



V. RESULTS 

A. Tube with quadratic cross section 

One of the few cases for which the analytical solution of the hydrodynamic problem Eqs. 
(@) " © i s known is Poiseuille flow, the flow through a linear tube with constant cross 
section. We will consider a tube with quadratic cross section, because here the geometry 
can be discretized on a cubic lattice without discretization error. 

We consider a tube directed along the t\ direction with quadratic cross section of side 
length B = 32a. We compare the velocity component ^1(22,^3) with its exactly known 
reference value v\ ei (x2, £3) given in Ref. | 40[| . Figure 0a shows the relative error {v\— v{ ei )/v[ ei 



for the LB solution with relaxation parameter r = 0.688. In Fig. we show the same 
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relative error for the LB solution with r = 1.0 (lower surface) and the FD solution (upper 
surface). The pressure gradient in the FD simulation was 2/33 while in the LB simulations 
it was around 10~ 4 . 

Around the center of the tube the analytical flow profile is very well recovered. Near the 
boundaries we find deviations which are biggest in the corners. The LB solution with r = 
0.688 underestimates the reference value t> ref while the solution with r = 1.0 overestimates 
the true value. Hence, the relaxation parameter r or equivalently the viscosity rj could 
be adjusted to find better agreement with the analytical velocity field. From Fig. ||] one 
expects to find a value 0.688 < r < 1.0 for which the numerical solution closely matches the 
analytical velocity profile. 

The computation time needed by the FD method which terminated when max^ e p | Av(x) — 
Vp(x)\ < 10~ 8 was 951 s on a DEC Alpha workstation. In LB simulations, the relative error 
of permeability was below 10~ 5 in 754 s on a Cray T3E for r = 1.0, but the simulations 
were continued for over 5000 s to make sure the saturation of the velocity fields. 

B. Cubic array of spheres 

To test the accuracy and efficiency of our two algorithms in a more complicated geometry 
with narrow constrictions we computed flow past a cubic array of spheres. This problem has 
become a reference system for checking hydrodynamic algorithms because accurate reference 
values for the permeability, and the drag coefficient, are available over a wide range of 
porosities jO], 0- 

The solution of this problem proceeds by solving the problem in a single unit cell of the 
cubic lattice. We generated six different unit cells of size L e {20a, 36a, 56a, 63a, 71a, 89a}. 
A sphere is placed at the center of each cell whose radius is chosen such that the porosity 
matches as close as possible to <fi = 0.15. Thus, the porosity is close to the porosity of the 
sandstones investigated later. 

The FD solution of the flow field was computed using periodic boundary conditions on 
those faces of the unit cell that are parallel to the macroscopic flow direction. On the faces 
of the unit cell perpendicular to the macroscopic flow we applied the conditions Eqs. (|l6|) - 
(^T|) with the standard pressure gradient p m = 1, p ont = —1- In the LB solution the flow field 
was computed using periodic boundary conditions in all directions. This difference of the 
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boundary conditions arises from the fact that in the LB simulations the density fluctuations 
around an average density are calculated while in the FD simulations the pressure field enters 
directly. 

Once the velocity field was known we calculated the permeability from Eq. (El). Following 



[IT | we then utilize the expression 

h 1 / T, \ 3 

(28) 



R 2 QttC d \R 

to obtain the reference value /c re f of the permeability from the drag coefficient Cd given in 



Ref . I^JJ . The radius R(4>, L) of the spheres depends on and L and is given implicitly by 



the expression 



8tt /R\ 3 n /R\ 2 vr 1 

T i - 3 ni • (29) 



To calculate the reference value fc re f we solve this equation and find (L/R) ~ 1.6011 for 
= 0.15. Using the drag coefficient C D = 1.020 x 10 3 @, we find k Tc{ /R 2 « 0.0002135. 



In Fig. |3] the relative error {k — /c re f)//c re f of the permeability is plotted as a function 
of the linear dimensionless system size L/a. With increasing resolution the results of both 
methods converge to each other and the error predominantly decreases. In the lattice- 
Boltzmann simulations the relaxation parameter r = 1.0 is seen to give regularly better 
results than r = 0.688. It thus appears that for r = 1.0 the effective location of the 
no-flow boundary is more satisfactory than for r = 0.688, and, consequently, the relative 
error is smaller for t — 1.0 even though the finite-size effects are similar in both cases. 
The deviation for L = 89a is negative and varies between 3% and 6% depending on the 
method. This discrepancy might possibly result mainly from discretization errors as there 
is a similar oscillatory trend in the LB as well as the FD results. Further work on larger 
systems is however needed to answer the question whether the discrepancy might also result 
from other sources. 

The curves in Fig. 3 are all nonmonotonous. This results most likely from the discretiza- 
tion of the cross-sectional area of the pore throats between the spheres. The curve of the 
discretized cross-sectional function of L/a shows a similar nonmonotonic behavior. 

Besides the accuracy, the demand of computation time is the second important charac- 
teristic of a numerical method. Comparison of the computation time of the LB and the FD 
algorithms is difficult. The FD method iterates the physical velocity field v and the physical 
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pressure field p. The iteration is terminated when v and p fulfill the Stokes equation Eq. 
(0) with a predefined accuracy. In the LB method on the other hand, the particle distribu- 
tions fi are iterated. The velocity and pressure fields are calculated from the final density 
distribution. In our implementation there is no deterministic stopping criterion, although 
such criterion could be included. In practice the calculations were terminated after a fixed 
number of iteration steps. This number was determined for each system by comparing runs 
of different length. 

All calculations with the FD code were performed on a Cray T3E-900/512 at the HLRS 
computing center of the University of Stuttgart with a peak performance of 461 Gflops. 
The LB code was run on a Cray T3E- 750/512 at Center for Scientific Computing (CSC) in 
Espoo, Finland, with a peak performance of 384 Gflops. In order to compare the run times 
for the two codes we took the actual time required to execute the program corrected by the 
ratio of peak performances. The run time required for the LB code was calculated from the 
number of iterations multiplied by a conversion factor. The conversion factor was 0.039448 
for L = 20a, 0.208096 for L = 36a, 0.628828 for L = 56a, 0.885556 for L = 63a, 1.275176 for 
L = 71a, and 2.21636 for L = 89a. It was determined by the wall time spent for one iteration 
step computed from averages over several 100-step iterations. The memory requirements of 
the two algorithms are different. The FD algorithm requires to store 8 numbers per lattice 
node in the version used here. The D3Q19 model used for the LB-algorithm requires to store 
19 numbers per lattice node 

In Fig. § we compare the time evolution of the numerical value of the permeability k 
for different system sizes. Plotted on the x axis is the total time in seconds needed on two 
Cray-processors. For L = 20a and L = 36a, both methods reach the final value of k in 
approximately the same time. For large L the LB method seems to be faster. Notice that 
the results of the LB simulations shown in Fig. ^ were performed for r = 0.688. For r = 1.0 
the simulations were about 45 % faster. 

The convergence of k(t) towards its asymptotic value is monotonic for the FD method, 
while in the LB case k(t) shows strong oscillations for all L. The reason for these oscillations 



is probably the slight compressibility error inherent in the model |21|, |22[. The horizontal 
line in Fig. ^ gives the reference value fc re f towards which the asymptotic values of both 
algorithms converge with L. 

We also compared the permeabilities given by the Navier-Stokes and Stokes (linearized) 
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versions of the LB method for the cubic arrays of spheres. Very much as expected, the two 
results were the same within the first seven digits (body force about 10~ 4 ). 

C. Three-dimensional sandstones 

We now apply both algorithms to the solution of the hydrodynamic problem Eqs. @ - 
(H) within the irregular pore space geometries of the whole experimental sample EX and the 
model samples DM, GF, and SA. 

The FD algorithm used the boundary conditions as described above. The iteration scheme 
was terminated when the condition max^ 6 p \ Av(x) — Vp(af)| < 10~ 6 for the dimensionless 
left hand side of Eq. (|2]) was fulfilled for the first time. Thus, the relative error e(ku)/kn 
of the diagonal elements of the permeability tensor is estimated to be smaller than 0.012 in 
the case of EX, and e(ku)/ku < 0.36 in the case of SA. The relative error for the samples 
DM and GF lies in between these two extreme values. 

In the LB simulations no flow boundary conditions were applied on the sample surfaces 
parallel to the main flow direction. At the inlet and outlet (i.e. the sample surfaces perpen- 
dicular to the main flow direction) an additional fluid layer with a thickness of 19 — 21 lattice 
spacings was added and then periodic boundary conditions were applied. The body force 
did not act in the additional fluid layer. These additional fluid layers increased the total 
number of lattice points by about 8 % in comparison with the FD method. The relaxation 
parameter was r = 0.688. The simulation stopped after a predefined number of iterations 
which was estimated to suffice for the permeability to converge. 

In Table |TJ we give the components of the permeability tensors for all four samples and 
for both algorithms. 

The permeability results confirm the predictions from a previous purely geometrical anal- 
ysis based on local porosity theory H]. The analysis in Ref. |7| emphasized the importance 
of local connectivity. The permeabilities are strongly correlated with the geometrical con- 
nectivity of the pore space, measured by means of the total fraction of percolating cells 
Ps(L). In accordance with our discussion of ps(60a) given in Table | of Section III, we find 
that the permeability of the original sandstone EX and that of the process model DM are in 
good agreement, while the permeabilities of the stochastic models GF and SA are an order 
of magnitude smaller. It seems as if the stochastic reconstruction models cannot reproduce 
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the high degree of geometrical connectivity present in the original sandstone. The recon- 
structed two-point correlation function lacks information about the geometrical connectivity 
of the pore space. A correct description of the geometrical connectivity is however an indis- 
pensable precondition for the correct dynamical connectivity that determines the transport 
properties. 

We now proceed to compare the numerically obtained value of the permeability °f 
the Fontainebleau sandstone EX with the experimental value k* — 1.3 D. Such a comparison 
requires a correction due to the difference between the porosity of the EX sample and the 
porosity 0* = 0.1484 of the original core sample on which the experiment was performed. 
There exists a well-known experimental correlation between porosity and permeability of 
Fontainebleau sandstone |IIJ]. This correlation is usually approximated in the form 



in which A and b are constants. In the porosity range of interest ~ 0.13. . .0.15, this 
correlation has b ~ 4, with however a large uncertainty due to the scatter in the measured 
results. Hence, we can extrapolate the numerically determined permeabilities k = (ku + 
^22 + ^33) /3 into the prediction 



where 0* is the previously defined porosity of the core sample and the constant A has dropped 
out. From Eq. (|31~1) we obtain k* = 1150 mD for the FD method and k* = 1015 mD for the 
LB algorithm. These values are surprisingly close to the experimental value k* = 1300 mD. 
Such an excellent agreement is not common. This will in fact be seen in the following when 
we determine the permeability of a subsample. 

We also checked by the LB methods the difference between the Navier-Stokes and Stokes 
permeability of sample EX. The relative difference was found to be 0.00036 for the parame- 
ters specified above, with an about 17% longer simulation time in the full Navier-Stokes case. 
The smallness of this difference only demonstrates that, for the small pressure differences 
considered here, we indeed are in the Stokes regime. 



k = A<p b , 



(30) 




(31) 
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D. Fine Graining 



We investigate the permeability of subsamples of the original samples for two reasons. 
Firstly this gives us an estimate of the magnitude of permeability fluctuations, and secondly 
it allows us to estimate the dependence of numerically obtained permeability on the lattice 
resolution a. The subsamples have dimensions of Mi x Mi x M3 = 100 x 100 x 100. To test 



the dependence of k on a, we apply a fine graining scheme P5| . The fine graining algorithm 
replaces each lattice point by n x n x n lattice points of the same phase with n G {2, 3, 4}. 
Thus, we get systems of dimensions Mi x M 2 x M 3 = n ■ 100 x n ■ 100 x n ■ 100 with lattice 
spacing a n = — /im. For each system the hydrodynamic problem Eqs. (0) - is solved, 
and the permeability k^ n ' is determined. 

Figure [5] shows k^ for a cubic subsample of EX with an applied pressure gradient in 
direction e%. The permeabilities obtained from the LB simulations are significantly higher 
for all r. This is due to the extra fluid layer with a thickness of 10 % of M\ outside the sample 
and the periodic boundary conditions used in the LB simulations. For this smaller piece of 
the sample, the effect of the boundary conditions is more significant than for the whole 
sample. We also performed with LB a couple of simulations with conditions similar to those 
used in FD. The pressure boundary condition was imitated using the body force combined 
with density and momentum averaging at the inlet and the outlet in the adjacent free fluid 
layer with thickness of one lattice spacing. No-flow boundary conditions were applied on 
the other cube sides. We found that the results of the LB and the FD method were again 
very close. This is indicated in Fig. ^| by two isolated points (filled square and star) at 
a\ = 7.5 /iin which were obtained from LB simulations with these boundary conditions for 
t = 0.688 and r = 1.0, respectively. However, a complete recomputation of the data of Fig. 
[5] would have exhausted the available computation time. Moreover, the results for the free 
fluid layer and for periodic boundary conditions allow us to estimate the influence of the 
boundary conditions on the result. 

For the original resolution the permeabilities obtained from the LB algorithm differ with 
varying r by nearly a factor 2. This again shows that the accuracy of the LB results for 
low-porosity (and low-discretization) systems strongly depends on the relaxation parameter. 
While for r = 0.6 the LB results are nearly independent of the lattice spacing, the changes in 
the permeability are drastically increased with increasing r. As we have already discussed, 
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in the permeability simulations, decrease in the relaxation parameter r (or viscosity) can be 
used to compensate an inadequate grid resolution (finite-size effects) within certain limits. 
By testing the structures of the subsamples for several grid resolutions and values of r, we 
chose r = 0.688 for the main simulations of the permeabilities. This is not necessarily the 
optimal choice and it differs from the results for pipe flow, but a more accurate calibration 
of r would require further simulations or other independent results. 

The permeability k^ for an infinite resolution of the lattice is obtained from a linear 
extrapolation of the data displayed in Fig. ||. Extrapolation of the FD results yields k^°°' = 
224 niD. Thus, we obtain an estimate for the relative error of our permeabilities in the case 
of EX, (jfeW - k^)/k^ « 0.33. For LB the extrapolated permeability is fc(°°) = 550 mD. 
Analogously, for r = 0.688 we obtain (jfeW - A; (oo) )/A; (1) « 0.15. For DM this error of LB is 
about the same size as for EX, and approximately 0.3 — 0.4 for SA and GF. 

The simulations of the subsamples show also that the magnitude of fluctuations in per- 
meability, within one and the same sample, can be 100 % even when the sample is extremely 
homogeneous. An effect of similar size may be induced by inaccurate boundary conditions. 
This is important when comparing two micropermeameter experiments in which different 
boundary conditions may have applied. 

When analyzing a subsample of sample SA, we encountered another difference between 
the FD and LB simulation. Geometrical analysis of the considered subsample of SA reveals 
that this subsample is not percolating in direction e\. We found nevertheless that its LB 
permeability is k^ = 50 mD for r = 0.688, and for the original resolution a\ = 7.5 /im. We 
attribute this result to diagonal leaks that are present in the LB model used. 

E. Velocity distributions 

Next we consider the velocity fields in more detail, and analyze the histograms of velocity 
magnitudes. Figure |6] shows the scaled distributions of the magnitude of the velocity \v(x)\ 
with f 6 P for samples EX and SA. The distributions were sampled using the solutions of 
the flow fields for an applied pressure gradient in direction e\. The higher permeability of 
the original sandstone EX is reflected by a higher probability density of regions with average 
flow velocity. The distribution of sample SA on the other hand exhibits a higher peak at 
\v\ =0, and it extends to higher velocities. The former observation indicates large stagnant 
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areas where no transport is taking place. The increased probability at high velocities may 
be related to a large number of narrow pore throats through which the fluid has to move. 
The velocity results for sample GF resemble closely to those of sample SA while those for 
sample DM resemble those for EX. 

Unlike the results for the permeabilities, which are close to each other, the velocity 
distributions reveal more significant differences between FD and LB calculations (see Fig. 
P). The distributions obtained from the LB solution exhibit a maximum at small velocities 
which is not present in the distributions obtained from the FD solution. The differences 
could perhaps be attributed to the slip velocities at the boundary due, e.g., to the bounce- 
back boundary condition and the diagonal leak flows in the LB algorithm, which could lead 
to a systematic deficit of zero velocities near the boundary. 

For SA the LB simulations show an additional spurious (double) peak at |#|/(|#|) ~ 1.5. 
Similar, although smaller peaks were also found for GF and EX. The locations of these peaks 
were found to be AP/2p{\v\), where AP is the body force. We attribute these peaks to the 
staggered momenta found in small closed pores, the number of which is very high for SA and 
GF. We have checked that these peaks disappear when time-averaged velocities are used, 
and a corresponding increase appears at zero velocity, which is the expected velocity for 
small closed pores. 

VI. CONCLUSION 

We have performed numerical micropermeametry on three dimensional porous mi- 
crostructures with a linear size of approximately 3 mm and a resolution of 7.5 /im. One 
of the samples has been a microtomographic image of Fontainebleau sandstone. Two of 
the samples were stochastic reconstructions with the same porosity, specific surface area, 
and two-point correlation function as the Fontainebleau sample. The fourth sample was a 
physical model, which mimics the processes of sedimentation, compaction and diagenesis of 
Fontainebleau sandstone. The permeabilities of these samples were determined by numeri- 
cally solving at low Reynolds numbers the appropriate Stokes equations in the pore spaces 
of the samples, using standard finite differences methods and the lattice-Boltzmann method. 
Our work shows that both methods, the LB method as well as standard FD methods, are 
applicable to the solution of the steady-state Stokes equation within the microstructure of a 



21 



three-dimensional porous medium. We investigated systems with sizes of up to 400 3 lattice 
points. The solution of even larger systems seems possible. Hence, numerical microperme- 
ametry is becoming a feasible technique for studying permeability fluctuations. 

An accurate, quantitative comparison of the two numerical methods is difficult due to the 
different approaches underlying these methods. The memory requirements of the algorithms 
used in this study differ by roughly a factor of 2.5. The FD algorithm requires to store 8 real 
numbers per lattice node, while the LB algorithm needs 19 real numbers per node. However, 



in this LB model 15 real numbers could also be used [19]. Considering the time consumption, 
both methods are quite similar. Our comparison has shown that there are some features in 
the standard application of the LB method, which need special attention. These include the 
r-dependence of the no-flow boundary, the compressibility of the fluid, staggered invariants, 
and diagonal leak flows. One should notice that these occasionally inconvenient features can 
usually be eliminated if so needed. The compressibility of the fluid can be eliminated for 



stationary flows |44|, and the effects of staggered invariants can be eliminated with proper 
averaging. The diagonal leak flow, which becomes noticeable for rough surfaces and near the 
percolation threshold and limits there the accuracy of the method, can also be eliminated 
by fairly straightforward means. The r-dependence of the no-flow boundaries can as well 
be eliminated, entirely by introducing a modified LB model [53], or almost entirely by using 
second-order boundaries as discussed above. We could of course have implemented here all 
these corrections in the LB code, and achieved thereby a more favourable comparison with 
the FD code and experiment, but we wanted to show the points of concern in a "standard" 
implementation of the LB method. 

Our results provide a quantitative comparison of various models for porous rocks. We 
show for the first time that stochastic reconstruction models for Fontainebleau are less 



accurate than originally believed f47|- In addition, our results for the permeabilities of 
the Fontainebleau sandstone and its models confirmed previous predictions 0] of a purely 
geometrical investigation of the same samples based on local porosity theory. The numerical 
value of k for the EX sample was found to be in good agreement with the experimental 
value. Nevertheless, one has to keep in mind that "numerically exact" results for the fluid 
permeability must be handled with care. We have shown here that they may depend on the 
numerical method, the boundary conditions, the size of the sample, and the resolution of the 
microstructure. Errors of as much as 100 % cannot typically been ruled out. In summary, 
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numerically exact determination of permeability is difficult to achieve. On the other hand, 
the same is true for the precision measurements in an experiment. 
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EX DM GF SA 



4> 


0.1355 


0.1356 


0.1421 


0.1354 


S V (M) [mm- 1 ] 


9.99 


10.30 


14.53 


11.04 


K V (M)[ mm.- 2 } 


-151 


-194 


-449 


-222 


T V (M) [mm" 3 ] 


-2159 


-2766 


4334 


14484 


fp [%] 


99.35 


99.23 


79.16 


62.73 


p 3 (60a) 


0.9561 


0.9647 


0.3255 


0.1695 



TABLE I: Geometrical characteristics of the four samples. Sy, Ky and Ty are specific surface, 
specific integral of mean curvature and specific integral of total curvature of the matrix phase, 
respectively, f p is the fraction of percolating pore lattice points, and ps(60a) is the probability of 
finding a cubic subblock of size L = 60a of the sample, which is percolating in all three directions. 
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FIG. 1: Spatial discretization of the velocity and pressure field on a MAC grid. The grey shaded 
cell lies in matrix phase, the white cell in pore space. 
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TABLE II: Permeability tensors of the Fontainebleau sandstone and its models. The values are 
given in mD. 
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FIG. 2: Relative error of the numerical solution of the velocity field for a Poiseuille flow through 
a tube with quadratic cross section of size B = 32a. The upper surface shows the FD solution, 
the lower surfaces the LB solutions with r = 1.0 and r = 0.688, respectively. As the FD solution, 
the LB solution with r = 1.0 overestimates the reference solution while the LB solution with 
r = 0.688 underestimates the references values. The reference values are calculated from the 
analytical solution given in Ref. ifjcj ]. 
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FIG. 3: Numerical results for the permeability k of a cubic array of spheres for different values of 
the lattice spacing a. The porosity is constant for all systems, <p = 0.15. The reference value k re { 
is taken from Ref. pi]] . 
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FIG. 4: Time evolution of the numerical solution of the permeability k for flow through a cubic 
array of spheres with porosity (f> = 0.15. The relaxation parameter was r = 0.688 for the LB 
simulations. 
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FIG. 5: Permeability k of a cubic subsample of the Fontainebleau sandstone (EX) for different 
values of the lattice spacing a and different relaxation parameters r. The size of the subsample 
is L = 750 //m. The LB simulation used an additional fluid layer and the periodic boundary 
conditions. The isolated data points with a = 7.5 fj,m are obtained from LB simulation without 
extra fluid layer, with momentum averaging in the inlet and outlet and with no-flow boundary 
conditions. The LB relaxations parameters were r = 1.0 and r = 0.688. 
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FIG. 6: Velocity distribution function -P(|i?|) of the Fontainebleau sandstone (EX) and the SA model 
sampled over the pore space. The distributions are scaled with the mean velocity (\v(x)\)^^. The 
inset shows a magnification of the distribution of EX for small velocities. 
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